home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / hlvector.lha / hl_vector / hjmin.cc < prev    next >
C/C++ Source or Header  |  1993-08-08  |  6KB  |  203 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *              Numerical Math Package
  6.  *        Hook-Jeevse multidimensional minimization
  7.  *
  8.  * Synopsis
  9.  *    double hjmin(b,h0,funct)
  10.  *    Vector& b             should be specified to contain the
  11.  *                    initial guess to the point of minimum.
  12.  *                    On return, it contains the location
  13.  *                    of the minimum the function has found
  14.  *    const Vector& h0        Initial values for steps along
  15.  *                    each direction
  16.  *    double f(const Vector& x)    Procedure to compute a function
  17.  *                    value at the specified point
  18.  *
  19.  * The function returns the value of the minimizing function f() at the
  20.  * point of the found minimum.
  21.  *
  22.  * Algorithm
  23.  *    Hook-Jeevse method of direct search for a function minimum
  24.  *    The method is of the 0. order (i.e. requiring no gradient computation)
  25.  *    See
  26.  *    B.Bondi. Methods of optimization. An Introduction - M.,
  27.  *    "Radio i sviaz", 1988 - 127 p. (in Russian)
  28.  *    
  29.  ************************************************************************
  30.  */
  31.  
  32.  
  33. #include "LinAlg.h"
  34. #include "math_num.h"
  35. #include <std.h>
  36.  
  37. #pragma implementation
  38.  
  39. /*
  40.  *------------------------------------------------------------------------
  41.  *        Class to operate on the points in the space (x,f(x))
  42.  */
  43.  
  44. class FPoint
  45. {
  46.   Vector& x;                // Point in the function domain
  47.   double fval;                // Function value at the point
  48.   double (*fproc)(const Vector& x);    // Procedure to compute the function
  49.                     // value
  50.   const short free_x_on_destructing;    // The flag is needed to free the 
  51.                     // dynamic memory allocated when
  52.                     // x rather than reference to x has
  53.                     // been given to construct FPoint
  54.  
  55. public:
  56.   FPoint(const Vector& b, double (*f)(const Vector& x));
  57.   FPoint(const FPoint& fp);
  58.   ~FPoint();
  59.  
  60.   FPoint& operator = (const FPoint& fp);
  61.  
  62.   double f() const            { return fval; }
  63.  
  64.   double examination(const Vector& h);    // Examine the function in the
  65.                     // neighborhood of the current point.
  66.                     // h defines the radius of the region
  67.  
  68.                     // Proceed in direction the function
  69.                     // seems decline
  70.   friend void update_in_direction(FPoint& from, FPoint& to);
  71.  
  72.                     // Decide whether the region embracing
  73.                     // the local min is small enough
  74.   int is_step_relatively_small(const Vector& h, const double tau);
  75. };
  76.  
  77.                 // Constructor FPoint from array b
  78. inline FPoint::FPoint(const Vector& b, double (*f)(const Vector& x))
  79.     : x(b), fproc(f), free_x_on_destructing(0)
  80. {
  81.   fval = (*fproc)(x);
  82. }
  83.  
  84.                 // Constructor by example
  85. FPoint::FPoint(const FPoint& fp)
  86.     : x(*(new Vector(fp.x))), free_x_on_destructing(1)
  87. {
  88.   x = fp.x;
  89.   fval = fp.fval;
  90.   fproc = fp.fproc;
  91. }
  92.                 // Destructor
  93. FPoint::~FPoint()
  94. {
  95.   if( free_x_on_destructing )
  96.     delete &x;
  97. }
  98.  
  99.                 // Assignment
  100. inline FPoint& FPoint::operator = (const FPoint& fp)
  101. {
  102.   x = fp.x;
  103.   fval = fp.fval;
  104.   fproc = fp.fproc;
  105.   return *this;
  106. }
  107.  
  108. /*
  109.  * Examine the function f in the vicinity of the current point
  110.  * by making tentative steps fro/back along each coordinate.
  111.  * Should function decrease, the point is updated to locate the
  112.  * new local min.
  113.  * Examination() returns the minimal function value found in
  114.  * the region.
  115.  *
  116.  */
  117. double FPoint::examination(const Vector& h)
  118. {
  119.   register int i;
  120.             // Perform a step along a coordinate
  121.   for(i=x.q_lwb(); i<=x.q_upb(); i++)
  122.   {
  123.     const double hi = h(i);
  124.     register double xi_old = x(i);      // Old value of x[i]
  125.     register double fnew;
  126.  
  127.     if( x(i) = xi_old + hi, fnew = (*fproc)(x), fnew < fval )
  128.       fval = fnew;            // Step caused f to decrease, OK
  129.     else if( x(i) = xi_old - hi, fnew = (*fproc)(x), fnew < fval )
  130.       fval = fnew;
  131.     else                // No function decline has been
  132.       x(i) = xi_old;                    // found along this coord, back up
  133.   }
  134.   return fval;
  135. }                                                
  136.  
  137.                 // Proceed in the direction the function
  138.                 // seems to fall
  139.                 // to_new = (to - from) + to
  140.                 // from = to (before modification)
  141. void update_in_direction(FPoint& from, FPoint& to)
  142. {
  143.   register int i;
  144.   for(i=(from.x).q_lwb(); i<=(from.x).q_upb(); i++)
  145.   {
  146.     register double t = (to.x)(i);
  147.     (to.x)(i)  += (t - (from.x)(i));
  148.     (from.x)(i) = t;
  149.   }
  150.   from.fval = to.fval;
  151.   to.fval = (*(to.fproc))(to.x);
  152. }
  153.  
  154.                 // Estimate whether point of minimum has
  155.                 // been located accurately enough
  156. int FPoint::is_step_relatively_small(const Vector& h, const double tau)
  157. {
  158.   register int it_is_small = 1;
  159.   register int i;
  160.   for(i=h.q_lwb(); it_is_small && i<=h.q_upb(); i++)
  161.     it_is_small &= ( h(i) /(1 + fabs(x(i))) < tau );
  162.   return it_is_small;
  163. }
  164.  
  165. /*
  166.  *------------------------------------------------------------------------
  167.  *                Root module
  168.  */
  169.  
  170. double hjmin(Vector& b, const Vector& h0, double (*ff)(const Vector& x))
  171. {
  172.                 // Function Parameters
  173.   const double tau = 10*EPSILON;        // Termination criterion
  174.   const double threshold = 1e-8;        // Threshold for the function
  175.                                         // decay to be treated as
  176.                                         // significant
  177.   const double step_reduce_factor = 10;   
  178.  
  179.   are_compatible(b,h0);
  180.  
  181.   FPoint pmin(b,ff);            // Point of min
  182.   FPoint pbase(pmin);            // Base point
  183.  
  184.   Vector h(h0); h = h0;            // Current steps along the axes
  185.  
  186.   for(;;)            // Main iteration loop
  187.   {                         // pmin is the approximation to min so far
  188.     if( pbase.examination(h) < pmin.f() - threshold )
  189.     {                         // Function value dropped significantly
  190.       do                                  // from pmin to the point pbase
  191.     update_in_direction(pmin,pbase);// Keep going in the same direction
  192.       while( pbase.examination(h) < pmin.f() - threshold ); // while it works
  193.       pbase = pmin;            // Save the best approx found
  194.     }
  195.     else                                   // Function didn't fall significantly
  196.       if(                               // upon wandering around pbas
  197.      h *= 1/step_reduce_factor,    // Try to reduce the step then
  198.      pbase.is_step_relatively_small(h,tau) )
  199.     return pmin.f();
  200.   }
  201. }
  202.  
  203.